home *** CD-ROM | disk | FTP | other *** search
/ Developer CD Series 2000 November: Tool Chest / Dev.CD Nov 00 TC Disk 1.toast / Sample Code / Devices and Hardware / Velocity Engine / VelEng Multiprecision / vfactor source / factor.c next >
Encoding:
C/C++ Source or Header  |  2000-09-28  |  15.4 KB  |  922 lines  |  [TEXT/CWIE]
  1. /**************************************************************
  2.  *
  3.  *    vfactor.c
  4.  *
  5.  *    AltiVec-based factoring program.
  6.  *  Created as extension of original factor.c project at 
  7.  *  Next Software, Inc.
  8.  *
  9.  *  This package is part of ongoing research in the
  10.  *    Advanced Computation Group, Apple Computer.
  11.  *    
  12.  *    c. 1999 Apple Computer, Inc.
  13.  *    All Rights Reserved.
  14.  *
  15.  *
  16.  *************************************************************/
  17.  
  18. /*
  19.     Disclaimer:    IMPORTANT:  This Apple software is supplied to you by Apple Computer, Inc.
  20.                 ("Apple") in consideration of your agreement to the following terms, and your
  21.                 use, installation, modification or redistribution of this Apple software
  22.                 constitutes acceptance of these terms.  If you do not agree with these terms,
  23.                 please do not use, install, modify or redistribute this Apple software.
  24.  
  25.                 In consideration of your agreement to abide by the following terms, and subject
  26.                 to these terms, Apple grants you a personal, non-exclusive license, under Apple’s
  27.                 copyrights in this original Apple software (the "Apple Software"), to use,
  28.                 reproduce, modify and redistribute the Apple Software, with or without
  29.                 modifications, in source and/or binary forms; provided that if you redistribute
  30.                 the Apple Software in its entirety and without modifications, you must retain
  31.                 this notice and the following text and disclaimers in all such redistributions of
  32.                 the Apple Software.  Neither the name, trademarks, service marks or logos of
  33.                 Apple Computer, Inc. may be used to endorse or promote products derived from the
  34.                 Apple Software without specific prior written permission from Apple.  Except as
  35.                 expressly stated in this notice, no other rights or licenses, express or implied,
  36.                 are granted by Apple herein, including but not limited to any patent rights that
  37.                 may be infringed by your derivative works or by other works in which the Apple
  38.                 Software may be incorporated.
  39.  
  40.                 The Apple Software is provided by Apple on an "AS IS" basis.  APPLE MAKES NO
  41.                 WARRANTIES, EXPRESS OR IMPLIED, INCLUDING WITHOUT LIMITATION THE IMPLIED
  42.                 WARRANTIES OF NON-INFRINGEMENT, MERCHANTABILITY AND FITNESS FOR A PARTICULAR
  43.                 PURPOSE, REGARDING THE APPLE SOFTWARE OR ITS USE AND OPERATION ALONE OR IN
  44.                 COMBINATION WITH YOUR PRODUCTS.
  45.  
  46.                 IN NO EVENT SHALL APPLE BE LIABLE FOR ANY SPECIAL, INDIRECT, INCIDENTAL OR
  47.                 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
  48.                 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  49.                 ARISING IN ANY WAY OUT OF THE USE, REPRODUCTION, MODIFICATION AND/OR DISTRIBUTION
  50.                 OF THE APPLE SOFTWARE, HOWEVER CAUSED AND WHETHER UNDER THEORY OF CONTRACT, TORT
  51.                 (INCLUDING NEGLIGENCE), STRICT LIABILITY OR OTHERWISE, EVEN IF APPLE HAS BEEN
  52.                 ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  53. */
  54.  
  55. /* include files */
  56.  
  57. #include <stdio.h>
  58. #include <math.h>
  59. #include <stdlib.h>
  60. #include <time.h>
  61. #include <events.h>
  62.  
  63. #define DO_RANDOM_CURVESEED    0 // use random curves (vs. iterate from start seed)
  64. #define VERBOSE    1              // print curve info and show progress
  65.  
  66. #ifdef _WIN32 
  67.  
  68. #include <process.h>
  69.  
  70. #endif
  71.  
  72. #include <string.h>
  73. #if __VEC__
  74. #include "vgiants.h"
  75. #else
  76. #include "origgiants.h"
  77. #endif
  78.  
  79.  
  80. /* definitions */
  81.  
  82. #define D 100
  83. #define NUM_PRIMES 6542 /* PrimePi[2^16]. */
  84.  
  85.  
  86. /* compiler options */
  87.  
  88. #ifdef _WIN32
  89. #pragma warning( disable : 4127 4706 ) /* disable conditional is constant warning */
  90. #endif
  91.  
  92.  
  93. /* global variables */
  94.  
  95. extern giant     scratch2;
  96. int             pr[NUM_PRIMES];
  97. giant             xr = NULL, xs = NULL, zs = NULL, zr = NULL, xorg = NULL,
  98.                 zorg = NULL, t1 = NULL, t2 = NULL, t3 = NULL, N = NULL,
  99.                 gg = NULL, An = NULL, Ad = NULL;
  100. giant             xb[D+2], zb[D+2], xzb[D+2];
  101. int             modmode = 0, Q, modcount = 0;
  102.  
  103.  
  104. /* function prototypes */
  105.  
  106. void        ell_even(giant x1, giant z1, giant x2, giant z2, giant An,
  107.                         giant Ad, giant N);
  108. void        ell_odd(giant x1, giant z1, giant x2, giant z2, giant exor,
  109.                         giant zor, giant N);
  110. void        ell_mul(giant xx, giant zz, int n, giant An, giant Ad, giant N);
  111. int         least_2(int n);
  112.  
  113. #if VERBOSE
  114. void        dot(void);
  115. #else
  116. #define dot()
  117. #endif
  118.  
  119. int            psi_rand();
  120.  
  121.  
  122. /**************************************************************
  123.  *
  124.  *    Functions
  125.  *
  126.  **************************************************************/
  127.  
  128.  
  129. int            
  130. psi_rand(
  131.     void
  132. )
  133. {
  134.     unsigned short    hi;
  135.     unsigned short    low;
  136.     time_t            tp;
  137.     int                result;
  138.  
  139.     time(&tp);
  140.     low = (unsigned short)rand();
  141.     hi = (unsigned short)rand();
  142.     result = ((hi << 16) | low) ^ ((int)tp);
  143.  
  144.     return (result & 0x7fffffff);
  145. }
  146.  
  147.  
  148. static void
  149. set_random_seed(
  150.     void
  151. )
  152. {
  153.     /* Start the random number generator at a new position. */
  154.     time_t        tp;
  155.  
  156.     time(&tp);
  157.     srand((int)tp + (int)TickCount());
  158. }
  159.  
  160.  
  161. static int
  162. isprime(
  163.     int     odd
  164. )
  165. {
  166.     int     j;
  167.     int     p;
  168.  
  169.     for (j=1; ; j++)
  170.     {
  171.         p = pr[j];
  172.         if (p*p > odd)
  173.             return(1);
  174.         if (odd % p == 0)
  175.             return(0);
  176.     }
  177. }
  178.  
  179.  
  180. static int
  181. primeq(
  182.     int                p
  183. )
  184. {
  185.     register int    j=3;
  186.  
  187.     if ((p&1)==0)
  188.         return ((p==2)?1:0);
  189.     if (j>=p)
  190.         return (1);
  191.     while ((p%j)!=0)
  192.     {
  193.         j+=2;
  194.         if (j*j>p)
  195.             return(1);
  196.     }
  197.     return(0);
  198. }
  199.  
  200.  
  201. static void
  202. s_modg(
  203.     giant        N,
  204.     giant        t
  205. )
  206. {
  207.     ++modcount;
  208.     switch (modmode)
  209.     {
  210.         case 0:
  211.             modg(N, t);
  212.             break;
  213.         case -1:
  214.             mersennemod(Q, t);
  215.             break;
  216.         case 1:
  217.             fermatmod(Q, t);
  218.             break;
  219.     }
  220. }
  221.  
  222.  
  223. static void
  224. reset_mod(
  225.     giant     x,
  226.     giant     N
  227. )
  228. /* Perform a divide (by the discovered factor) and switch back
  229.    to non-Fermat-non-Mersenne (i.e. normal) mod mode. */
  230. {
  231.     divg(x, N);
  232.     modmode = 0;
  233. }
  234.  
  235. void
  236. ell_even(
  237.     giant     x1,
  238.     giant     z1,
  239.     giant     x2,
  240.     giant     z2,
  241.     giant     An,
  242.     giant     Ad,
  243.     giant     N
  244. )
  245. {
  246.     gtog(x1, t1);
  247.     addg(z1, t1);
  248.     squareg(t1);
  249.     s_modg(N, t1);
  250.     gtog(x1, t2);
  251.     subg(z1, t2);
  252.     squareg(t2);
  253.     s_modg(N, t2);
  254.     gtog(t1, t3);
  255.     subg(t2, t3);
  256.     gtog(t2, x2);
  257.     mulg(t1, x2);
  258.     gshiftleft(2, x2);
  259.     s_modg(N, x2);
  260.     mulg(Ad, x2);
  261.     s_modg(N, x2);
  262.     mulg(Ad, t2);
  263.     gshiftleft(2, t2);
  264.     s_modg(N, t2);
  265.     gtog(t3, t1);
  266.     mulg(An, t1);
  267.     s_modg(N, t1);
  268.     addg(t1, t2);
  269.     mulg(t3, t2);
  270.     s_modg(N, t2);
  271.     gtog(t2,z2);
  272. }
  273.  
  274.  
  275. void
  276. ell_odd(
  277.     giant     x1,
  278.     giant     z1,
  279.     giant     x2,
  280.     giant     z2,
  281.     giant     exor,
  282.     giant     zor,
  283.     giant     N
  284. )
  285. {
  286.     gtog(x1, t1);
  287.     subg(z1, t1);
  288.     gtog(x2, t2);
  289.     addg(z2, t2);
  290.     mulg(t1, t2);
  291.     s_modg(N, t2);
  292.     gtog(x1, t1);
  293.     addg(z1, t1);
  294.     gtog(x2, t3);
  295.     subg(z2, t3);
  296.     mulg(t3, t1);
  297.     s_modg(N, t1);
  298.     gtog(t2, x2);
  299.     addg(t1, x2);
  300.     squareg(x2);
  301.     s_modg(N, x2);
  302.     gtog(t2, z2);
  303.     subg(t1, z2);
  304.     squareg(z2);
  305.     s_modg(N, z2);
  306.     mulg(zor, x2);
  307.     s_modg(N, x2);
  308.     mulg(exor, z2);
  309.     s_modg(N, z2);
  310. }
  311.  
  312.  
  313. void
  314. ell_mul(
  315.     giant             xx,
  316.     giant             zz,
  317.     int             n,
  318.     giant             An,
  319.     giant             Ad,
  320.     giant             N
  321. )
  322. {
  323.     unsigned int     c = (unsigned int)0x80000000;
  324.  
  325.     if (n==1)
  326.         return;
  327.     if (n==2)
  328.     {
  329.         ell_even(xx, zz, xx, zz, An, Ad, N);
  330.         return;
  331.     }
  332.     gtog(xx, xorg);
  333.     gtog(zz, zorg);
  334.     ell_even(xx, zz, xs, zs, An, Ad, N);
  335.  
  336.     while((c&n) == 0)
  337.     {
  338.         c >>= 1;
  339.     }
  340.  
  341.     c>>=1;
  342.     do
  343.     {
  344.         if (c&n)
  345.         {
  346.             ell_odd(xs, zs, xx, zz, xorg, zorg, N);
  347.             ell_even(xs, zs, xs, zs, An, Ad, N);
  348.         }
  349.         else
  350.         {
  351.             ell_odd(xx, zz, xs, zs, xorg, zorg, N);
  352.             ell_even(xx, zz, xx, zz, An, Ad, N);
  353.         }
  354.         c >>= 1;
  355.     } while(c);
  356. }
  357.  
  358.  
  359.  
  360. /* From R. P. Brent, priv. comm. 1996:
  361. Let s > 5 be a pseudo-random seed (called $\sigma$ in the Tech. Report),
  362.  
  363.     u/v = (s^2 - 5)/(4s)
  364.  
  365. Then starting point is (x_1, y_1) where
  366.  
  367.     x_1 = (u/v)^3
  368. and
  369.     a = (v-u)^3(3u+v)/(4u^3 v) - 2
  370. */
  371.  
  372. static void
  373. choose12(
  374.     giant     x,
  375.     giant     z,
  376.     int     k,
  377.     giant     An,
  378.     giant     Ad,
  379.     giant     N
  380. )
  381. {
  382.     itog(k, zs);
  383.     gtog(zs, xs);
  384.     squareg(xs);
  385.     itog(5, t2);
  386.     subg(t2, xs);
  387.     s_modg(N, xs);
  388.     addg(zs, zs);
  389.     addg(zs, zs);
  390.     s_modg(N, zs);
  391.     gtog(xs, x);
  392.     squareg(x);
  393.     s_modg(N, x);
  394.     mulg(xs, x);
  395.     s_modg(N, x);
  396.     gtog(zs, z);
  397.     squareg(z);
  398.     s_modg(N, z);
  399.     mulg(zs, z);
  400.     s_modg(N, z);
  401.  
  402.     /* Now for A. */
  403.     gtog(zs, t2);
  404.     subg(xs, t2);
  405.     gtog(t2, t3);
  406.     squareg(t2);
  407.     s_modg(N, t2);
  408.     mulg(t3, t2);
  409.     s_modg(N, t2);  /* (v-u)^3. */
  410.     gtog(xs, t3);
  411.     addg(t3, t3);
  412.     addg(xs, t3);
  413.     addg(zs, t3);
  414.     s_modg(N, t3);
  415.     mulg(t3, t2);
  416.     s_modg(N, t2);  /* (v-u)^3 (3u+v). */
  417.     gtog(zs, t3);
  418.     mulg(xs, t3);
  419.     s_modg(N, t3);
  420.     squareg(xs);
  421.     s_modg(N, xs);
  422.     mulg(xs, t3);
  423.     s_modg(N, t3);
  424.     addg(t3, t3);
  425.     addg(t3, t3);
  426.     s_modg(N, t3);
  427.     gtog(t3, Ad);
  428.     gtog(t2, An);  /* An/Ad is now A + 2. */
  429. }
  430.  
  431.  
  432. static void
  433. ensure(
  434.     int    q
  435. )
  436. {
  437.     int     nsh, j;
  438.  
  439.     N = newgiant(0);
  440.     if(!q)
  441.     {
  442.         gread(N,stdin);
  443.         q = bitlen(N) + 1;
  444.     }
  445.     nsh = q/4; /* Allowing (easily) enough space per giant,
  446.                     since N is about 2^q, which is q bits, or
  447.                    q/16 shorts.  But squaring, etc. is allowed,
  448.                     so we need at least q/8, and we choose q/4
  449.                     to be conservative. */
  450.     if (!xr)
  451.         xr = newgiant(nsh);
  452.     if (!zr)
  453.         zr = newgiant(nsh);
  454.     if (!xs)
  455.         xs = newgiant(nsh);
  456.     if (!zs)
  457.         zs = newgiant(nsh);
  458.     if (!xorg)
  459.         xorg = newgiant(nsh);
  460.     if (!zorg)
  461.         zorg = newgiant(nsh);
  462.     if (!t1)
  463.         t1 = newgiant(nsh);
  464.     if (!t2)
  465.         t2 = newgiant(nsh);
  466.     if (!t3)
  467.         t3 = newgiant(nsh);
  468.     if (!gg)
  469.         gg = newgiant(nsh);
  470.     if (!An)
  471.         An = newgiant(nsh);
  472.     if (!Ad)
  473.         Ad = newgiant(nsh);
  474.     for (j=0;j<D+2;j++)
  475.     {
  476.         xb[j] = newgiant(nsh);
  477.         zb[j] = newgiant(nsh);
  478.         xzb[j] = newgiant(nsh);
  479.     }
  480. }
  481.  
  482. static int
  483. bigprimeq(
  484.     giant     x
  485. )
  486. {
  487.     itog(1, t1);
  488.     gtog(x, t2);
  489.     subg(t1, t2);
  490.     itog(5, t1);
  491.     
  492.     powermodg(t1, t2, x);
  493.     if (isone(t1))
  494.         return(1);
  495.     return(0);
  496. }
  497.  
  498. #if VERBOSE
  499. void
  500. dot(
  501.     void
  502. )
  503. {
  504.     printf(".");
  505.     fflush(stdout);
  506. }
  507. #endif
  508. /**************************************************************
  509.  *
  510.  *    Main Function
  511.  *
  512.  **************************************************************/
  513.  
  514. int main(
  515.     int    argc,
  516.     char     *argv[]
  517. )
  518. {
  519.     int     j, k, C, nshorts, cnt, count,
  520.             limitbits = 0, pass, npr, rem;
  521.     long    B;
  522.     int     randmode = 0;
  523.     long     cntseed;
  524.  
  525.  
  526.     if (!strcmp(argv[argc-1], "-r"))
  527.     {
  528.         randmode = 1;
  529.         if (argc > 4)
  530.           /* This segment only takes effect in random mode. */
  531.             limitbits = atoi(argv[argc-2]);
  532.     }
  533.     else
  534.     {
  535.         randmode = 0;
  536.     }
  537.  
  538.     modmode = 0;
  539.     if (argc > 2)
  540.     {
  541.         modmode = atoi(argv[1]);
  542.         Q = atoi(argv[2]);
  543.     }
  544.     if (modmode==0)
  545.         Q = 0;
  546.     ensure(Q);
  547.     if (modmode)
  548.     {
  549.         itog(1, N);
  550.         gshiftleft(Q, N);
  551.         itog(modmode, t1);
  552.         addg(t1, N);
  553.     }
  554.     pr[0] = 2;
  555.     for (k=0, npr=1;; k++)
  556.     {
  557.         if (primeq(3+2*k))
  558.         {
  559.             pr[npr++] = 3+2*k;
  560.             if (npr >= NUM_PRIMES)
  561.                 break;
  562.         }
  563.     }
  564.  
  565.  
  566.     if (randmode == 0)
  567.     {
  568.         printf("Sieving...\n");
  569.         fflush(stdout);
  570.         for (j=0; j < NUM_PRIMES; j++)
  571.         {
  572.             gtog(N, t1);
  573.             rem = idivg(pr[j], t1);
  574.             if (rem == 0)
  575.             {
  576.                 printf("%d ", pr[j]);
  577.                 gtog(t1, N);
  578.                 if (isone(N))
  579.                 {
  580.                     printf("\n");
  581.                     goto exit;
  582.                 }
  583.                 else
  584.                 {
  585.                     printf("* ");
  586.                     fflush(stdout);
  587.                 }
  588.                 --j;
  589.             }
  590.         }
  591.  
  592.  
  593.         if (bigprimeq(N))
  594.         {
  595.             gout(N);
  596.             goto exit;
  597.         }
  598.  
  599.         printf("\n");
  600.         printf("Commencing Pollard rho...\n");
  601.         fflush(stdout);
  602.         itog(1, gg);
  603.         itog(3, t1); itog(3, t2);
  604.  
  605.         for (j=0; j < 15000; j++)
  606.         {
  607.             if((j%100) == 0)
  608.             {
  609.                 dot();
  610.                 gcdg(N, gg);
  611.                 if (!isone(gg))
  612.                     break;
  613.             }
  614.             squareg(t1);
  615.             iaddg(2, t1);
  616.             s_modg(N, t1);
  617.             squareg(t2);
  618.             iaddg(2, t2);
  619.             s_modg(N, t2);
  620.             squareg(t2);
  621.             iaddg(2, t2);
  622.             s_modg(N, t2);
  623.             gtog(t2, t3);
  624.             subg(t1, t3);
  625.             absg(t3);
  626.             mulg(t3, gg);
  627.             s_modg(N, gg);
  628.         }
  629.         gcdg(N, gg);
  630.  
  631.         if ((gcompg(N,gg) != 0) && (!isone(gg)))
  632.         {
  633.             fprintf(stdout,"\n");
  634.             gout(gg);
  635.             reset_mod(gg, N);
  636.             if (isone(N))
  637.             {
  638.                 printf("\n");
  639.                 goto exit;
  640.             }
  641.             else
  642.             {
  643.                 printf("* ");
  644.                 fflush(stdout);
  645.             }
  646.             if (bigprimeq(N))
  647.             {
  648.                 gout(N);
  649.                 goto exit;
  650.             }
  651.         }
  652.  
  653.         printf("\n");
  654.         printf("Commencing Pollard (p-1)...\n");
  655.         fflush(stdout);
  656.         itog(1, gg);
  657.         itog(3, t1);
  658.         for (j=0; j< NUM_PRIMES; j++)
  659.         {
  660.             cnt = (int)(8*log(2.0)/log(1.0*pr[j]));
  661.             if (cnt < 2)
  662.                 cnt = 1;
  663.             for (k=0; k< cnt; k++)
  664.             {
  665.                 powermod(t1, pr[j], N);
  666.             }
  667.             itog(1, t2);
  668.             subg(t1, t2);
  669.             mulg(t2, gg);
  670.             s_modg(N, gg);
  671.  
  672.             if (j % 100 == 0)
  673.             {
  674.                 dot();
  675.                 gcdg(N, gg);
  676.                 if    (!isone(gg))
  677.                     break;
  678.             }
  679.         }
  680.         gcdg(N, gg);
  681.         if ((gcompg(N,gg) != 0) && (!isone(gg)))
  682.         {
  683.             fprintf(stdout,"\n");
  684.             gout(gg);
  685.             reset_mod(gg, N);
  686.             if (isone(N))
  687.             {
  688.                 printf("\n");
  689.                 goto exit;
  690.             }
  691.             else
  692.             {
  693.                 printf("* ");
  694.                 fflush(stdout);
  695.             }
  696.             if (bigprimeq(N))
  697.             {
  698.                 gout(N);
  699.                 goto exit;
  700.             }
  701.         }
  702.     } /* This is the end of (randmode == 0) */
  703.  
  704.     printf("\n");
  705.     printf("Commencing ECM...\n");
  706.     fflush(stdout);
  707.  
  708.     if (randmode)
  709.         set_random_seed();
  710.     pass = 0;
  711.     
  712.     
  713. #if !DO_RANDOM_CURVESEED
  714.     cntseed = 1438858426;        // initial curve seed point
  715. #endif
  716.  
  717.     
  718.     while (++pass)
  719.     {
  720.         if (randmode == 0)
  721.         {
  722.             if (pass <= 3)
  723.             {
  724.                 B = 1000;
  725.             }
  726.             else if (pass <= 10)
  727.             {
  728.                 B = 10000;
  729.             }
  730.             else if (pass <= 100)
  731.             {
  732.                 B = 100000L;
  733.             } else
  734.             {
  735.                 B = 1000000L;
  736.             }
  737.         }
  738.         else
  739.         {
  740.             B = 2000000L;
  741.         }
  742.         C = 50*((int)B);
  743.  
  744.         /* Next, choose curve with order divisible by 16 and choose
  745.          *    a point (xr/zr) on said curve.
  746.          */
  747.  
  748.         /* Order-div-12 case. 
  749.          * cnt = 8020345;   Brent's parameter for stage one discovery
  750.          * of 27-digit factor of F_13.
  751.          */
  752. #if DO_RANDOM_CURVESEED
  753.         cnt = psi_rand(); /* cnt = 8020345; */
  754. #else
  755.         cnt = cntseed++;
  756. #endif
  757.         choose12(xr, zr, cnt, An, Ad, N);
  758. #if VERBOSE         
  759.         printf("Choosing curve %d, with s = %d, B = %d, C = %d:\n", pass,cnt, B, C);   fflush(stdout);
  760. #endif        
  761.         cnt = 0;
  762.         nshorts = 1;
  763.         count = 0;
  764.         for (j=0;j<nshorts;j++)
  765.         {
  766.             ell_mul(xr, zr, 1<<16, An, Ad, N);
  767.             ell_mul(xr, zr, 3*3*3*3*3*3*3*3*3*3*3, An, Ad, N);
  768.             ell_mul(xr, zr, 5*5*5*5*5*5*5, An, Ad, N);
  769.             ell_mul(xr, zr, 7*7*7*7*7*7, An, Ad, N);
  770.             ell_mul(xr, zr, 11*11*11*11, An, Ad, N);
  771.             ell_mul(xr, zr, 13*13*13*13, An, Ad, N);
  772.             ell_mul(xr, zr, 17*17*17, An, Ad, N);
  773.         }
  774.         k = 19;
  775.         while (k<B)
  776.         {
  777.             if (isprime(k))
  778.             {
  779.                 ell_mul(xr, zr, k, An, Ad, N);
  780.                 if (k<100)
  781.                     ell_mul(xr, zr, k, An, Ad, N);
  782.                 if (cnt++ %100==0)
  783.                     dot();
  784.             }
  785.             k += 2;
  786.         }
  787.         count = 0;
  788.  
  789.         gtog(zr, gg);
  790.         gcdg(N, gg);
  791.         if ((!isone(gg))&&(bitlen(gg)>limitbits))
  792.         {
  793.             fprintf(stdout,"\n");
  794.             gwriteln(gg, stdout);
  795.             fflush(stdout);
  796.             reset_mod(gg, N);
  797.             if (isone(N))
  798.             {
  799.                 printf("\n");
  800.                 goto exit;
  801.             }
  802.             else
  803.             {
  804.                 printf("* ");
  805.                 fflush(stdout);
  806.             }
  807.             if (bigprimeq(N))
  808.             {
  809.                  gout(N);
  810.                  goto exit;
  811.             }
  812.             continue;
  813.         }
  814.         else
  815.         {
  816. #if VERBOSE
  817.             printf("\n");            
  818.             fflush(stdout);
  819. #endif            
  820.         }
  821.  
  822.         /* Continue;  Invoke, to test Stage 1 only. */
  823.         k = ((int)B)/D;
  824.         gtog(xr, xb[0]);
  825.         gtog(zr, zb[0]);
  826.         ell_mul(xb[0], zb[0], k*D+1 , An, Ad, N);
  827.         gtog(xr, xb[D+1]);
  828.         gtog(zr, zb[D+1]);
  829.         ell_mul(xb[D+1], zb[D+1], (k+2)*D+1 , An, Ad, N);
  830.  
  831.         for (j=1; j <= D; j++)
  832.         {
  833.             gtog(xr, xb[j]);
  834.             gtog(zr, zb[j]);
  835.             ell_mul(xb[j], zb[j], 2*j , An, Ad, N);
  836.             gtog(zb[j], xzb[j]);
  837.             mulg(xb[j], xzb[j]);
  838.             s_modg(N, xzb[j]);
  839.         }
  840.         modcount = 0;
  841. #if VERBOSE
  842.         printf("\nCommencing second stage, curve %d...\n",pass); fflush(stdout);
  843. #endif        
  844.         count = 0;
  845.         itog(1, gg);
  846.  
  847.         while (1)
  848.         {
  849.             gtog(zb[0], xzb[0]);
  850.             mulg(xb[0], xzb[0]);
  851.             s_modg(N, xzb[0]);
  852.             mulg(zb[0], gg);
  853.             s_modg(N,gg); /* Accumulate. */
  854.             for (j = 1; j < D; j++)
  855.             {
  856.                 if (!isprime(k*D+1+ 2*j))
  857.                     continue;
  858.  
  859.                 /* Next, accumulate (xa - xb)(za + zb) - xa za + xb zb. */
  860.                 gtog(xb[0], t1);
  861.                 subg(xb[j], t1);
  862.                 gtog(zb[0], t2);
  863.                 addg(zb[j], t2);
  864.                 mulg(t1, t2);
  865.                 s_modg(N, t1);
  866.                 subg(xzb[0], t2);
  867.                 addg(xzb[j], t2);
  868.                 s_modg(N, t2);
  869.                 --modcount;
  870.                 mulg(t2, gg);
  871.                 s_modg(N, gg);
  872.                 if((++count)%1000==0)
  873.                     dot();
  874.             }
  875.  
  876.             k += 2;
  877.             if(k*D > C)
  878.                 break;
  879.             gtog(xb[D+1], xs);
  880.             gtog(zb[D+1], zs);
  881.             ell_odd(xb[D], zb[D], xb[D+1], zb[D+1], xb[0], zb[0], N);
  882.             gtog(xs, xb[0]);
  883.             gtog(zs, zb[0]);
  884.         }
  885.  
  886.         gcdg(N, gg);
  887.         if((!isone(gg))&&(bitlen(gg)>limitbits))
  888.         {
  889.             fprintf(stdout,"\n");
  890.             gwriteln(gg, stdout);
  891.             fflush(stdout);
  892.             reset_mod(gg, N);
  893.             if (isone(N))
  894.             {
  895.                 printf("\n");
  896.                 goto exit;
  897.             }
  898.             else
  899.             {
  900.                 printf("* ");
  901.                 fflush(stdout);
  902.             }
  903.             if (bigprimeq(N))
  904.             {
  905.                 gout(N);
  906.                 goto exit;
  907.             }
  908.             continue;
  909.         }
  910.  
  911. #if VERBOSE
  912.         printf("\n");
  913.         fflush(stdout);
  914. #endif        
  915.     }
  916.  
  917. exit:
  918.  
  919.     return 0;
  920. }
  921.  
  922.